<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of funcrs</title>
  <meta name="keywords" content="funcrs">
  <meta name="description" content="Cross approximation of a function of a TT-tensor, Method 1">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../index.html">Home</a> &gt;  <a href="#">tt2</a> &gt; <a href="index.html">@tt_tensor</a> &gt; funcrs.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../index.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for tt2/@tt_tensor&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>funcrs
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>Cross approximation of a function of a TT-tensor, Method 1</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>function [y]=funcrs(tt,fun,eps,y,nswp,varargin) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">Cross approximation of a function of a TT-tensor, Method 1
   [Y]=FUNCRS(TT,FUN,EPS,Y,NSWP)
   Computes approximation to the function FUN(TT) with accuracy EPS
   Auxiliary parameters:  Y (initial approximation), NSWP 
   (number of sweeps in the method 
   Much faster then usual cross by vectorized computation of subtensors

   
 TT-Toolbox 2.2, 2009-2012

This is TT Toolbox, written by Ivan Oseledets et al.
Institute of Numerical Mathematics, Moscow, Russia
webpage: http://spring.inm.ras.ru/osel

For all questions, bugs and suggestions please mail
ivan.oseledets@gmail.com
---------------------------</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="../../tt2/@qtt_tucker/diag.html" class="code" title="function [qt]=diag(qt)">diag</a>	Diagonal of a matrix or diagonal matrix from a vector in QTT-Tucker</li><li><a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>	Frobenius norm of the QTT-Tucker</li><li><a href="../../tt2/@tt_matrix/conj.html" class="code" title="function [b]=conj(a)">conj</a>	Complex conjugate of a TT-matrix</li><li><a href="../../tt2/@tt_matrix/core.html" class="code" title="function [tt] = core(tt1,varargin)">core</a>	Converts TT-matrix to TT1 cell-array format</li><li><a href="../../tt2/@tt_matrix/diag.html" class="code" title="function [tt]=diag(tm)">diag</a>	Extract the diagonal of the TT-matrix</li><li><a href="../../tt2/@tt_matrix/norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>	Matrix norm of the TT-matrix</li><li><a href="../../tt2/@tt_matrix/size.html" class="code" title="function [sz] = size(tt)">size</a>	Mode sizes of the TT-matrix</li><li><a href="conj.html" class="code" title="function [tt1]=conj(tt)">conj</a>	Compute a complex conjugate of TT-tensor</li><li><a href="core.html" class="code" title="function [tt] = core(tt1,varargin)">core</a>	Converts TT-tensor TT1 to old-cell array format.</li><li><a href="diag.html" class="code" title="function [tm]=diag(tt)">diag</a>	Constructs diagonal TT-matrix from TT-tensor</li><li><a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>	Frobenius norm of the TT-tensor</li><li><a href="qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>	Left and right orthogonalization of the TT-format</li><li><a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>	Reshape of the TT-tensor</li><li><a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>	Mode sizes of the TT-tensor</li><li><a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>	Maximal volume submatrix in an tall matrix</li><li><a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>	Truncation by absolution precision in Frobenius norm</li><li><a href="../../tt2/core/reort.html" class="code" title="function [y]=reort(u,uadd)">reort</a>	Golub-Kahan reorthogonalization</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
</ul>
<!-- crossreference -->



<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [y]=funcrs(tt,fun,eps,y,nswp,varargin)</a>
0002 <span class="comment">%Cross approximation of a function of a TT-tensor, Method 1</span>
0003 <span class="comment">%   [Y]=FUNCRS(TT,FUN,EPS,Y,NSWP)</span>
0004 <span class="comment">%   Computes approximation to the function FUN(TT) with accuracy EPS</span>
0005 <span class="comment">%   Auxiliary parameters:  Y (initial approximation), NSWP</span>
0006 <span class="comment">%   (number of sweeps in the method</span>
0007 <span class="comment">%   Much faster then usual cross by vectorized computation of subtensors</span>
0008 <span class="comment">%</span>
0009 <span class="comment">%</span>
0010 <span class="comment">% TT-Toolbox 2.2, 2009-2012</span>
0011 <span class="comment">%</span>
0012 <span class="comment">%This is TT Toolbox, written by Ivan Oseledets et al.</span>
0013 <span class="comment">%Institute of Numerical Mathematics, Moscow, Russia</span>
0014 <span class="comment">%webpage: http://spring.inm.ras.ru/osel</span>
0015 <span class="comment">%</span>
0016 <span class="comment">%For all questions, bugs and suggestions please mail</span>
0017 <span class="comment">%ivan.oseledets@gmail.com</span>
0018 <span class="comment">%---------------------------</span>
0019 
0020 
0021 <span class="comment">%PARAMETERS SECTION</span>
0022 rmin=1; 
0023 dropsweeps = 4;
0024 kick_rank=10;
0025 verb=true;
0026 <span class="keyword">if</span> (~isempty(y))
0027    yold=y;
0028 <span class="keyword">end</span>
0029 <span class="comment">%For this procedure we need only local (!) indices, since it</span>
0030 <span class="comment">%is more or less equivalent to orthogonalization;</span>
0031 d=tt.d; 
0032 ps=tt.ps;
0033 <a href="core.html" class="code" title="function [tt] = core(tt1,varargin)">core</a>=tt.core;
0034 n=tt.n;
0035 r=tt.r;
0036 
0037 ry=y.r;
0038 psy=y.ps;
0039 cry=y.core;
0040 phi=cell(d+1,1);
0041 phi{d+1}=1;
0042 phi{1}=1;
0043 <span class="comment">%Warmup procedure: orthogonalize from right-to-left &amp; maxvol</span>
0044 pos1=psy(d);
0045 <span class="comment">%Cores i</span>
0046 <span class="keyword">for</span> i=d-1:-1:1
0047    cr=cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1);
0048    cr2=cry(pos1-ry(i)*n(i)*ry(i+1):pos1-1);
0049    cr2=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[ry(i)*n(i),ry(i+1)]);
0050    cr=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr,[ry(i+1),n(i+1)*ry(i+2)]);
0051    cr=cr.';
0052    [cr,rm]=<a href="qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(cr,0); 
0053    ry(i+1)=<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(cr,2);
0054    ind=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(cr); 
0055    r1=cr(ind,:);
0056    cr=cr/r1;
0057    cr=cr.';
0058    cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=cr(:);
0059    pos1=pos1-ry(i)*n(i)*ry(i+1);
0060    cr2=cr2*(r1*rm).'; 
0061    cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=cr2(:);
0062    <span class="comment">%Take phi matrix; convolve from right with current cors of V</span>
0063    cr0=<a href="core.html" class="code" title="function [tt] = core(tt1,varargin)">core</a>(ps(i+1):ps(i+2)-1);
0064    cr0=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr0,[r(i+1)*n(i+1),r(i+2)]); 
0065    cr0=cr0*phi{i+2}; <span class="comment">%cr0 is now r(i)*n(i)*ry(i+1);</span>
0066    cr0=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr0,[r(i+1),n(i+1)*ry(i+2)]);
0067    phi{i+1}=cr0(:,ind); 
0068  <span class="comment">%  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);</span>
0069 <span class="comment">%y.core=cry;</span>
0070 <span class="comment">%y.r=ry;</span>
0071 <span class="comment">%y.ps=psy;</span>
0072 <span class="comment">%keyboard</span>
0073 
0074    <span class="comment">%Orthogonalization &amp; maxvol is &quot;local operation&quot; (touches two</span>
0075    <span class="comment">%cores)</span>
0076    <span class="comment">%Computation of elements requries storage of phi matrices</span>
0077 <span class="keyword">end</span>
0078 <span class="comment">%Truncate cry</span>
0079 <span class="comment">%pos1=pos1-n(1)*r(2);</span>
0080 cry=cry(pos1:numel(cry));
0081 y.core=cry;
0082 y.r=ry;
0083 y.ps=psy;
0084 <span class="comment">%keyboard</span>
0085 swp=1;
0086 
0087 yold=[];
0088 not_converged = true;
0089 pos1=1;
0090 last_sweep = false;
0091 <span class="keyword">while</span> ( swp &lt; nswp &amp;&amp; not_converged )
0092     max_er=0;
0093 cry_old=cry;
0094 psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0095 pos1=1;
0096  <span class="keyword">for</span> i=1:d-1
0097      <span class="comment">%fprintf('i=%d \n',i);</span>
0098      <span class="comment">%We care for two cores, with number i &amp; number i+1, and use</span>
0099      <span class="comment">%psi(i) and psi(i+2) as a basis; also we will need to recompute</span>
0100      <span class="comment">%psi(i+1)</span>
0101      ps1=phi{i}; ps2=phi{i+2};
0102      <span class="comment">%Take current core; convolve it with</span>
0103      <span class="comment">%core=core(ps</span>
0104      <span class="comment">%Compute (!) superblock and function (!) of it</span>
0105      cr1=<a href="core.html" class="code" title="function [tt] = core(tt1,varargin)">core</a>(ps(i):ps(i+1)-1);
0106      cr2=<a href="core.html" class="code" title="function [tt] = core(tt1,varargin)">core</a>(ps(i+1):ps(i+2)-1);
0107      cr1=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr1,[r(i),n(i)*r(i+1)]);
0108      cr1=ps1*cr1;
0109      cr2=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[r(i+1)*n(i+1),r(i+2)]);
0110      cr2=cr2*ps2;
0111      cr1=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr1,[ry(i)*n(i),r(i+1)]);
0112      cr2=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[r(i+1),n(i+1)*ry(i+2)]);
0113      cr=cr1*cr2;
0114      cr=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
0115      cr=fun(cr); <span class="comment">%Elements are evaluated here!</span>
0116      cr=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
0117      <span class="comment">%Check for local approximation of cr for the error</span>
0118      cry1=cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1);
0119 
0120      cry2=cry_old(psy(i+1):psy(i+2)-1);
0121      
0122      cry1=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cry1,[ry(i)*n(i),ry(i+1)]);
0123      cry2=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cry2,[ry(i+1),n(i+1)*ry(i+2)]);
0124      appr=cry1*cry2;
0125      er=<a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(appr-cr,<span class="string">'fro'</span>)/<a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(cr,<span class="string">'fro'</span>);
0126      max_er=max(er,max_er);
0127      <span class="comment">%Compute SVD of cr</span>
0128 
0129      <span class="keyword">if</span> (mod(swp,dropsweeps)~=0)&amp;&amp;(swp&gt;1)&amp;&amp;(~last_sweep)
0130          [u,s,v]=svd(cr-appr,<span class="string">'econ'</span>);
0131      <span class="keyword">else</span>
0132          [u,s,v]=svd(cr,<span class="string">'econ'</span>);
0133      <span class="keyword">end</span>;
0134      s=<a href="diag.html" class="code" title="function [tm]=diag(tt)">diag</a>(s);
0135      r2=<a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>(s,eps*<a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(cr,<span class="string">'fro'</span>)/sqrt(d-1));
0136      s=s(1:r2); u=u(:,1:r2); v=v(:,1:r2);
0137      
0138      v=<a href="conj.html" class="code" title="function [tt1]=conj(tt)">conj</a>(v)*<a href="diag.html" class="code" title="function [tm]=diag(tt)">diag</a>(s);
0139 
0140      <span class="keyword">if</span> (mod(swp,dropsweeps)~=0)&amp;&amp;(swp&gt;1)&amp;&amp;(~last_sweep)
0141          u = [cry1, u];
0142          v = [cry2.', v];
0143          [u,rv]=<a href="qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(u,0);
0144          v = v*(rv.');         
0145      <span class="keyword">else</span>
0146          <span class="keyword">if</span> (~last_sweep)
0147              <span class="comment">%Kick rank of u</span>
0148              ur=randn(<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(u,1),kick_rank);
0149              <span class="comment">%Orthogonalize ur to u by Golub-Kahan reorth</span>
0150              u=<a href="../../tt2/core/reort.html" class="code" title="function [y]=reort(u,uadd)">reort</a>(u,ur);
0151              radd=<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(u,2)-r2;
0152              <span class="keyword">if</span> ( radd &gt; 0 )
0153                  vr=zeros(<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(v,1),radd);
0154                  v=[v,vr];
0155              <span class="keyword">end</span>
0156          <span class="keyword">end</span>;
0157      <span class="keyword">end</span>;
0158      
0159      <span class="comment">%vr=randn(size(v,1),kick_rank);</span>
0160      <span class="comment">%v=reort(v,vr);</span>
0161      <span class="comment">%radd=size(v,2)-rnew;</span>
0162      <span class="comment">%if ( radd &gt; 0 )</span>
0163      <span class="comment">%    ur=zeros(size(u,1),radd);</span>
0164      <span class="comment">%    u=[u,ur];</span>
0165      <span class="comment">%end</span>
0166      <span class="comment">%rnew=rnew+radd;</span>
0167 
0168      
0169      
0170      
0171      <span class="comment">%u=[u,uadd]; v=[v,vadd];</span>
0172      <span class="comment">%[u,ru]=qr(u,0);</span>
0173      <span class="comment">%v=v*(ru.');</span>
0174      r2=<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(u,2);
0175      
0176      ind=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(u);
0177      r1=u(ind,:); 
0178      u=u/r1; v=v*r1'; v=v.';
0179      <span class="comment">%Compute new phi</span>
0180      ry(i+1)=r2;
0181      cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);
0182      pos1=pos1+ry(i)*n(i)*ry(i+1);
0183      cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:);
0184      
0185      <span class="comment">%Now we have to: cr1 with phi from the left</span>
0186      phi{i+1}=cr1(ind,:); <span class="comment">%phi{i+1}=phi{i+1};</span>
0187      
0188      
0189  <span class="keyword">end</span>
0190   <span class="comment">%Truncate local memory</span>
0191   cry=cry(1:pos1+ry(d)*n(d)*ry(d+1)-1);
0192   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0193 y.core=cry;
0194 y.r=ry;
0195 y.ps=psy;
0196 <span class="comment">%keyboard;</span>
0197 cry_old=cry;
0198 <span class="comment">%m=numel(cry);</span>
0199 <span class="comment">%cry=[zeros(size(cry)),cry];</span>
0200 <span class="comment">%pos1=pos1+m;</span>
0201 <span class="comment">%cry2=cry_old(psy(i+1):psy(i+2)-1);</span>
0202 cry=cry(psy(d):psy(d+1)-1); <span class="comment">%Start--only two cores</span>
0203  <span class="keyword">for</span> i=d-1:-1:1
0204      <span class="comment">%fprintf('i=%d \n',i);</span>
0205      <span class="comment">%We care for two cores, with number i &amp; number i+1, and use</span>
0206      <span class="comment">%psi(i) and psi(i+2) as a basis; also we will need to recompute</span>
0207      <span class="comment">%psi(i+1)</span>
0208      ps1=phi{i}; ps2=phi{i+2};
0209      <span class="comment">%Take current core; convolve it with</span>
0210      <span class="comment">%core=core(ps</span>
0211      <span class="comment">%Compute (!) superblock and function (!) of it</span>
0212      cr1=<a href="core.html" class="code" title="function [tt] = core(tt1,varargin)">core</a>(ps(i):ps(i+1)-1);
0213      cr2=<a href="core.html" class="code" title="function [tt] = core(tt1,varargin)">core</a>(ps(i+1):ps(i+2)-1);
0214      cr1=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr1,[r(i),n(i)*r(i+1)]);
0215      cr1=ps1*cr1;
0216      cr2=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[r(i+1)*n(i+1),r(i+2)]);
0217      cr2=cr2*ps2;
0218      cr1=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr1,[ry(i)*n(i),r(i+1)]);
0219      cr2=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[r(i+1),n(i+1)*ry(i+2)]);
0220      cr=cr1*cr2;
0221      cr=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
0222      cr=fun(cr); <span class="comment">%Elements are evaluated here!</span>
0223      cr=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
0224      <span class="comment">%Check for local approximation of cr for the error</span>
0225      <span class="comment">%cry1=cry(pos1-ry(i)*n(i)*ry(i+1):pos1-1);</span>
0226      cry1=cry_old(psy(i):psy(i+1)-1);
0227      <span class="comment">%cry2=cry(ry(:ry(i+1)*n(i+1)*ry(i+2));</span>
0228      <span class="comment">%cry1=cry(1:ry(i)*n(i)*ry(i+1));</span>
0229      <span class="comment">%cry2=cry(ry(i)*n(i)*ry(i+1)+1:ry(i)*n(i)*ry(i+1)+ry(i+1)*n(i+1)*ry(i+2));</span>
0230      cry2=cry(1:ry(i+1)*n(i+1)*ry(i+2));
0231      cry(1:ry(i+1)*n(i+1)*ry(i+2))=[];
0232      <span class="comment">%cry2=cry_old(psy(i+1):psy(i+2)-1);</span>
0233      <span class="comment">%cry(1:(ry(i+1)*n(i+1)*ry(i+2)))=[]; %Delete both of first cores</span>
0234      <span class="comment">%cry(1:ry(i)*n(i)*ry(i+1)+ry(i+1)*n(i+1)*ry(i+2))=[];</span>
0235      cry1=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cry1,[ry(i)*n(i),ry(i+1)]);
0236      cry2=<a href="reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cry2,[ry(i+1),n(i+1)*ry(i+2)]);
0237      appr=cry1*cry2;
0238      er=<a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(appr-cr,<span class="string">'fro'</span>)/<a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(cr,<span class="string">'fro'</span>);
0239      <span class="comment">%er</span>
0240      max_er=max(er,max_er);
0241           
0242      <span class="keyword">if</span> (mod(swp,dropsweeps)~=0)&amp;&amp;(swp&gt;1)&amp;&amp;(~last_sweep)
0243          [u,s,v]=svd(cr-appr,<span class="string">'econ'</span>);
0244      <span class="keyword">else</span>
0245          <span class="comment">%Compute SVD of cr</span>
0246          [u,s,v]=svd(cr,<span class="string">'econ'</span>);
0247      <span class="keyword">end</span>;
0248      s=<a href="diag.html" class="code" title="function [tm]=diag(tt)">diag</a>(s);
0249      r2=<a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>(s,eps*<a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(cr,<span class="string">'fro'</span>)/sqrt(d-1));
0250      s=s(1:r2); u=u(:,1:r2); v=<a href="conj.html" class="code" title="function [tt1]=conj(tt)">conj</a>(v(:,1:r2));
0251      <span class="comment">%Make it standard</span>
0252      u=u*<a href="diag.html" class="code" title="function [tm]=diag(tt)">diag</a>(s);
0253      
0254      
0255      <span class="keyword">if</span> (mod(swp,dropsweeps)~=0)&amp;&amp;(swp&gt;1)&amp;&amp;(~last_sweep)
0256          u = [cry1, u];
0257          v = [cry2.', v];
0258          [v,rv]=<a href="qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(v,0);
0259          u = u*(rv.');
0260      <span class="keyword">else</span>
0261          <span class="keyword">if</span> (~last_sweep)
0262              <span class="comment">%Kick rank</span>
0263              vr=randn(<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(v,1),kick_rank);
0264              v=<a href="../../tt2/core/reort.html" class="code" title="function [y]=reort(u,uadd)">reort</a>(v,vr);
0265              radd=<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(v,2)-r2;
0266              <span class="keyword">if</span> ( radd &gt; 0 )
0267                  ur=zeros(<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(u,1),radd);
0268                  u=[u,ur];
0269              <span class="keyword">end</span>
0270          <span class="keyword">end</span>;
0271      <span class="keyword">end</span>;
0272 
0273      
0274      
0275      r2=<a href="size.html" class="code" title="function [sz] = size(tt,dim)">size</a>(v,2);
0276      ind=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(v);
0277      r1=v(ind,:); 
0278      v=v/r1; v=v.';
0279      u=u*r1';
0280      <span class="comment">%Compute new phi;</span>
0281      ry(i+1)=r2;
0282      u=u(:); u=u'; v=v(:); v=v';
0283      cry=[u,v,cry];
0284      <span class="comment">%keyboard;</span>
0285      <span class="comment">%We need new memory for</span>
0286      <span class="comment">%cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:); %Here memory has to be (?) enlarged</span>
0287      <span class="comment">%if ( pos1 &lt;= ry(i)*n(i)*ry(i+1) ) %Have to enlarge memory</span>
0288      <span class="comment">%  cry=cry(pos1:numel(cry));</span>
0289      <span class="comment">%  cry=[u(:)',cry];</span>
0290      <span class="comment">%  pos1=ry(i)*n(i)*ry(i+1)+1;</span>
0291      <span class="comment">%else</span>
0292      <span class="comment">%   pos1=pos1-ry(i)*n(i)*ry(i+1);</span>
0293       <span class="comment">%  cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);</span>
0294      <span class="comment">%end</span>
0295      <span class="comment">%Now we have to: cr1 with phi from the left</span>
0296      phi{i+1}=cr2(:,ind); phi{i+1}=phi{i+1};
0297  
0298  <span class="keyword">end</span>
0299  <span class="comment">%Truncate local memory</span>
0300  <span class="comment">%cry=cry(pos1:numel(cry));</span>
0301 
0302 <span class="comment">% keyboard;</span>
0303   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0304 
0305 y.core=cry;
0306 y.r=ry;
0307 y.ps=psy;
0308 <span class="comment">%keyboard;</span>
0309 <span class="keyword">if</span> ( isempty(yold) )
0310  yold=y;
0311  er_nrm=1;
0312 <span class="keyword">else</span>
0313    
0314    er_nrm=<a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(yold-y)/<a href="norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(y);
0315    yold=y;
0316 <span class="keyword">end</span>
0317 <span class="keyword">if</span> ( verb )
0318  fprintf(<span class="string">'sweep=%d, er=%3.2e er_nrm=%3.2e \n'</span>,swp,max_er,er_nrm);
0319 <span class="keyword">end</span>
0320 <span class="keyword">if</span> (last_sweep)
0321     <span class="keyword">break</span>;
0322 <span class="keyword">end</span>;
0323 <span class="keyword">if</span> (er_nrm&lt;eps)
0324     last_sweep=true;
0325 <span class="keyword">end</span>;
0326 swp=swp+1;
0327 <span class="keyword">end</span>     
0328   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0329 
0330 y.core=cry;
0331 y.r=ry;
0332 y.ps=psy;
0333 
0334 <span class="keyword">end</span></pre></div>
<hr><address>Generated on Wed 08-Feb-2012 18:20:24 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>